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ABSTRACT 

A  FINITE  ELEMENT  ANALYSIS  OF  A 

SINGLE  LAP  SHEAR  ADHESIVE  JOINT 

WITH  AND  WITHOUT  AN  EDGE  CRACK 

by 

THOMAS  HOLM  CHRISTENSEN 


Submitted  to  the  Department  of  Materials  Science  and  Engineering  and  to  the 
Department  of  Ocean  Engineering  on  January  21,  1976  in  partial  fulfillment  of  the 
requirements  for  the  degrees  of  Master  of  Science  in  Materials  Science  and  Engineering,  and 
Master  of  Science  in  Ocean  Engineering. 


The  single  lap  shear  adhesive  joint  is  analyzed  both  with  and  without  a  short  edge  crack 
using  the  finite  element  method.  The  method  employed  is  more  accurate  than  previous 
numerical  solutions  due  to  the  use  of  a  more  accurate  element  (assumed  stress  hybrid  element 
and  crack  tip  superelement)  and  a  finer  mesh  than  previously  applied.  The  results  of  the 
analysis  without  a  crack  shows  a  large  variation  in  stress  across  the  adhesive  not  previously 
reported.  The  results  also  show  previous  determinations  of  stress  at  the  end  of  lap  to  be 
conservative.  The  results  also  specifically  address  longitudinal  stress  in  the  adherend  which 
have  been  discussed  experimentally,  but  not  published  numerically.  The  analysis  of  the  single 
lap  joint  with  a  short  edge  crack  represents  the  first  application  of  the  general  nine  node 
singular  hybrid  stress  "super  element"  to  a  real  problem.  Results  indicate  a  power  lav/  relation 
beLween  stress  intensity  factors  kj  and  kn,  and  the  modulus  ratio.  It  is  shown  that  kn  is  not 
only  of  the  same  magnitude  as  kJt  but  actually  slightly  larger,  which  indicates  that  a  failure 
mechanism  must  be  combined  mode.  The  analysis  indicates  no  difference  between  the 
uncracked  joint  and  the  cracked  one  (crack  length  2.5  adhesive  thicknesses)  at  distances  greater 
than  ten  adhesive  thicknesses  away  from  the  crack  tip. 
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Section  I:  Introduction 

This  research  was  undertaken  to  investigate  the  stress  distribution  in  a  single  lap  shear 
adhesive  bonded  joint.  The  analysis  of  the  problem  was  done  using  the  finite  element  method 
(FEM).  The  intent  was  to  determine  the  stress  field,  in  the  configuration  shown  in  Figure  1, 
both  with  no  crack  present  and  with  a  single  short  edge  crack.  In  particular,  this  research 
investigates  the  variation  in  stress  across  the  adhesive  layer  without  a  crack,  and  then 
analyzes  the  changes  in  the  solution  when  a  short  edge  crack  is  introduced.  The  short  crack 
(2.5  adhesive  thicknesses)  was  chosen  so  that  the  constraint  conditions  which  were  found  to  be 
stable  (i.e.,  no  kinematic  modes)  for  the  uncracked  case  would  still  be  appropriate  for  the  joint 
with  a  short  crack.  The  crack  length  used  herein  idealizes  a  shrinkage  or  residual  stress 
induced  edge  crack. 

To  fully  investigate  the  lap  shear  adhesive  bonded  joint  would  require  a  vast  parametric 
study  which  systematically  investigates  the  effect  of  geometric  variables,(e.g.,  overlap  length, 
adhesive  thickness,  adherend  thickness,  adherend  shape  in  the  overlap  region,  and  the  adhesive 
shape  at  the  end  of  lap,  crack  length,  and  location  of  crack(s)),  as  well  as  material  properties, 
(e.g.,  Young's  modulus  ratio  for  isotropic  adherend  and  adhesive,  different  Poisson  ratios  for 
adherend  and  adhesive  and  anisotropic  materials.)  This  research  concentrates  on  the  effect  of 
the  Young's  modulus  ratio  over  a  range  of  one  to  one  thousand  with  other  variables  held 
constant  and  over  varied  adhesive  thicknesses  with  other  variables  held  constant. 

Section  II  is  a  brief  critical  review  of  the  literature  on  adhesive  bonded  joints  as  it  relates 
to  this  research.  Section  III  is  a  brief  discussion  of  the  FEM.  Section  IV  is  a  discussion  of  the 
singular  crack  tip  element  or  super  element  used  in  the  research  reported  herein  in  the  analysis 


of  the  joints  with  short  edge  cracks.  Section  V  is  a  presentation  of  the  results  and  discussiun 
of  this  research  on  uncracked  adhesive  bonded  lap  shear  joints.  Section  VI  is  a  presentation  of 
the  results  for  an  adhesive  bonded  lap  shear  joint  with  a  single  edge  crack  and  a  comparison 
with  the  uncracked  case  as  presented  in  Section  V.  Section  VII  is  a  summary  of  the  results  of 
the  previous  sections  and  a  presentation  of  the  conclusions  to  be  drawn  from  this  research. 
Section  VIII  contains  suggestions  for  further  work  to  be  undertaken  on  lap  shear  joints,  using 
the  FEM. 

Section  II:  Literature  Review 

The  literature  on  adhesive  bonded  joints  can  be  divided  into  three  broad  categories. 
First,  there  are  the  "classical"  solutions  which,  more  or  less,  attempt  a  closed  form  solution  after 
making  simplifying  assumptions.  Second,  there  are  the  numerical  techniques  of  finite 
differences  and  the  FEM.  Third,  there  are  the  experimental  investigations  of  joint  strength, 
fracture  mechanisms,  and  rough  stress  distribution  approximation. 

i.  Classical  Solutions 

Goland  and  Reissner  [1]  were  the  first  to  include  the  effects  of  bending  with  those  of 
differential  straining  in  their  solution  to  the  adhesive  bonded  lap  joint  in  tension.  Their 
solution  assumes  the  adhesive  layer  to  be  very  thin  relative  to  other  typical  dimensions  in  the 
problem  and,  therefore,  ignores  this  thickness  in  calculating  the  stresses.  The  solution  assumes 
linear  normal  stress  and  parabolic  shear  distribution,  thereby  reducing  the  problem  to  that  of  a 
cylindrically  bent  plate  or  thin  beam  theory.   The  primary  weakness  in  their  solution  is  that  it 
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is  assumed  that  the  stresses  do  not  vary  across  the  adhesive  thickness.  The  error  in  assuming  a 
constant  state  of  stress  will  be  shown  in  Section  V,  when  the  result  of  this  research  is  presented. 
The  Goland  and  Reissner  approximations,  however,  have  been  widely  used  as  a  measure  of  the 
validity  of  experimental  studies,  such  as,  by  the  U.S.  Forest  Service  [2].  Although  modifications 
have  been  made  to  the  Goland  and  Reissner  stress  equations,  the  shape  of  the  stress  curves  are 
similar.  For  this  reason,  the  Goland  and  Reissner  equations  are  included  as  Appendix  I,  and  a 
comparison  of  their  solution  with  the  current  research  is  presented  for  a  typical  geometry  in 
Figure  2  (tearing  stress)  and  Figure  3  (shear  stress). 

Cornell  [3]  used  Goland  and  Reissner's  work  as  a  basis  and,  his  solution  has  the  same 
error  in  that  he  also  assumes  a  constant  state  of  stress  across  the  adhesive  layer.  Cornell  treats 
the  adherends  as  independent  beams  connected  by  an  infinite  number  of  shear  and  tension 
springs.  The  research  reported  herein  finds  that  Cornell's  solution  to  the  resulting  differential 
equations  greatly  underestimates  the  longitudinal  stress  in  the  adherend.  High  longitudinal 
stress  is  due  to  the  localized  bending  effect  near  the  end  of  the  lap.  This  has  been  shown  to 
exist  experimentally  by  Wang,  et  al.  [4],  and  others,  and  to  be  sufficiently  high  to  cause 
plastic  flow  and  even  failure  in  the  adherend. 

Not  all  of  the  "classical"  solutions  are  as  old  as  Goland  and  Reissner's  work  and  Cornell's 
work.  Schijve  [5]  has  done  a  thin  beam  analysis  of  lap  shear  problems  to  determine  the  effect 
of  clamping  conditions  and  permanent  deformation  on  secondary  bending.  Schijve,  as  in 
previously  reported  classical  works,  assumes  a  thin  adhesive  layer  with  a  constant  stress  state 
across  its  thickness.  His  work  indicates  that  permanent  set,  or  local  plastic  deformation,  has  a 
greater  effect  than  clamping  conditions  on  secondary  bending  (i.e.,  stress  at  the  end  of  lap). 
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This  again  points  to  the  importance  of  determining  the  longitudinal  stress  in  the  adherend, 
however,  very  few  papers  discuss  longitudinal  stresses  at  all. 

Chang  and  Muki  [6]  have  solved  the  lap  shear  problem  by  reducing  the  problem  to  a 
pair  of  Fredholm  integral  equations  of  the  second  kind,  which  they  then  solve  for  bond  line 
stresses  and  stress  concentration  factors  at  the  end  of  the  bonded  region.  The  solution  is 
dependent  on  the  angle  of  rotation  of  the  joint,  which  is  determined  by  the  Goland  and 
Reissner  solution.  Chang  and  Muki  also  assume  an  infinitessimally  thin  adhesive  layer  and 
traction  and  displacement  vectors  continuous  across  the  bond.  As  would  be  expected  with  the 
same  basic  assumptions  as  Goland  and  Reissner,  the  results  are  similar  including  a  constant 
state  of  stress  across  the  bond. 

Srinivas  [7]  has  divided  the  lap  into  regions  where  the  adhesive  layer  is  treated  as  an 
infinite  number  of  elastic  tension  and  shear  springs.  This  is  similar  to  the  method  reported  by 
Cornell.  Srinivas,  however,  does  not  use  thin  beam  theory.  Instead,  he  expands  the 
longitudinal  displacement  (u)  and  the  transverse  displacement  (w)  up  to  the  second  order  in  the 
transverse  coordinate  (z),  where  the  undetermined  coefficients  are  functions  of  the  longitudinal 
coordinate  (x)  only.  These  assumed  displacements  are  then  used  with  the  governing  equations 
of  elasticity  with  appropriate  boundary  conditions  within  each  region  to  obtain  a  solution  for 
strains,  which,  in  turn,  determines  stresses  through  the  stress-strain  relations.  The  formulation 
is  sufficiently  involved  that  a  digital  computer  is  required  to  solve  the  resulting  linear 
equations.  Because  the  displacement  functions  are  limited  to  the  second  order,  this  solution  is 
no  better  than  a  finite  element  formulation  with  a  single  hybrid  element  (discussed  in  Section 
III)  or  a  single  eight  node  displacement  model  quadrilateral  element  (found  in  introductory 
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texts  on  FEM  [8,9]),  across  the  adherend  thickness  with  similar  boundary  conditions  at  the 
bond.  Either  of  the  latter  would  be  considered  unreasonably  coarse,  yet  both  are  of  the  same 
degree  of  accuracy  as  that  reported  by  Srinivas.  There  are  no  "classical"  treatments  of  the  lap 
shear  adhesive  bonded  joint  with  a  crack.  This  is  not  surprising  since  the  classical  solutions 
assume  infinitessimal  adhesive  thickness  and,  therefore,  cannot  address  a  crack  in  a  region 
which  is  not  directly  treated  in  the  formulation  of  the  problem. 

ii.  Numerical  Techniques 

Wooley  and  Carver  [10]  used  the  FEM  to  analyze  the  adhesive  bonded  lap  shear  joint. 
Their  work  has  several  basic  flaws.  First,  the  physical  geometry  of  the  adhesive  bonded  lap 
joint,  (i.e.,  width  much  greater  than  overlap  length  or  thickness)  makes  plane  strain,  as  used  by 
all  previously  mentioned  investigators,  appropriate,  yet  Wooley  and  Carver  assume  plane  stress. 
Second,  they  use  quadrilateral  elements  comprised  of  four  triangular  subelements  which 
produce  constant  stress  within  each  element,  and  less  accurate  results.  This  shortcoming  could 
be  reduced  by  a  more  imaginative  use  of  mesh  arrangement.  Wooley  and  Carver  use  equal 
mesh  spacing  across  the  thickness  and  along  the  lap  region  which  causes  them  to  completely 
miss  the  very  high  stress  concentration  that  exists  at  the  end  of  lap  for  low  adherend  to 
adhesive  modulus  ratios  due  to  the  geometric  and  material  discontinuities. 

Harrison  and  Harrison  [11]  also  apply  the  FEM  to  the  lap  shear  joint,  and  likewise,  their 
research  has  several  flaws.  First,  and  foremost,  they  assume  a  symmetry  about  the  midplane  of 
the  bond  which  does  not  physically  exist.  They  hypothesize  that  the  real  lap  shear  problem  can 
be  approximated  by  a  linear  combination  of  uniform  displacement  boundary  conditions  on  the 
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bond  line.  None  of  the  previous  literature  leads  one  to  believe  that  this  is  a  reasonable 
assumption.  Second,  their  formulation  neglects  rotation  and,  therefore,  could  never  predict 
incipient  failure  in  the  longitudinal  direction  due  to  adherend  yielding.  Third,  and  of  less 
importance,  they  use  all  constant  strain  elements.  This  requires  far  more  degrees  of  freedom 
and,  therefore,  computation  time,  than  other  elements  discussed  in  Section  III. 

Pirvics  [12]  investigates  the  lap  shear  problem  using  a  finite  difference  technique.  The 
physical  dimensions  used  in  Pirvics'  analysis  appear  to  unduly  restrict  rotation  of  the  system. 
That  is  to  say,  a  fixed  end  condition  is  imposed  on  one  adherend  at  a  distance  of  only  one  inch 
from  the  joint.  Also,  an  extremely  short  overlap  length  of  0.12  inch  is  used,  well  out  of  the 
normal  commercially  used  range  of  0.50-1.50  inches.  The  finite  difference  technique  is 
inherently  a  slowly  convergent  method.  The  finite  difference  method  also  has  great  difficulty 
handling  boundaries.  A  boundary  requires  a  transformation  of  the  actual  shape  so  that 
method  of  images  can  approximate  the  real  problems  discontinuous  shape  at  the  boundary. 
Because  of  the  concentration  of  stresses  near  a  boundary  in  the  physical  problem,  and  the 
difficulty  in  treating  boundary  regions  by  the  finite  differences  method  noted  above,  the 
validity  of  using  a  finite  differences  technique  is  in  question.  Pirvics  did  not  specifically  state 
how  the  boundaries  are  to  be  treated. 

Adams  and  Peppiatt  [13]  use  the  FEM  in  analyzing  both  single  and  double  lap  shear 
adhesive  bonded  joints.  The  only  single  lap  shear  configuration  for  which  they  present  results 
has  a  square  ended  adhesive  layer,  which  does  not  closely  approximate  a  real  single  lap  joint. 
In  their  analysis  of  double  lap  joints,  they  allow  a  fillet  "spew"  which  is  more  realistic.  That  is 
to  say,  a  real  adhesive  joint  when  put  under  some  pressure  at  time  of  bonding  will  tend  to 
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force  some  adhesive  out  of  the  lap  region  and  form  a  fillet.  Their  analysis  uses  all  constant 
strain  triangular  elements  which  have  shortcomings  as  noted  above.  They  do  not  address 
longitudinal  stresses  in  the  adherend,  nor  do  they  address  a  variation  in  stress  across  the 
adhesive  layer. 

Although  some  literature  does  exist  on  applications  of  the  FEM  to  the  adhesive  bonded 
joint  containing  a  crack,  this  problem  and  a  discussion  of  this  literature  will  be  postponed  to 
Section  IV.  None  of  the  existing  literature  has  used  the  singular  hybrid  stress  element  in  the 
analysis  of  the  lap  shear  problem  containing  a  crack  in  the  joint. 

iii.  Experimental  Techniques 

Experimental  studies  of  lap  joint  configuration  face  the  difficult  task  of  trying  to 
determine  the  stress  distribution  in  a  layer  only  five  to  ten  thousandths  of  an  inch  thick.  For 
this  reason,  most  experimental  studies  use  large  models  of  the  joint,  and  they  apply  the 
following  modeling  assumptions  [14]: 

a)  Maintain  the  ratio  of  adherend  thickness  to  adhesive  thickness  of 
the  actual  joint, 

b)  Maintain  the  same  ratio  of  adherend  shear  modulus  to  adhesive 
shear  modulus,  or  the  general  ratio  of  mechanical  properties  of 
adherend  to  adhesive,  as  in  the  actual  joint, 

c)  Be  wide  enough  to  ensure  a  condition  of  plane  strain  in  the  plane 
of  the  adhesive, 

d)  Simulate  a  true  adhesive  bond  between  the  adherend  and  adhesive. 
It  is  important  to  simulate  the  discontinuous  nature  of  the 
mechanical  properties  at  the  adherend-adhesive  interface. 

The  experimental  literature  on  adhesive  bonded  lap  joints  can  be  broken   into  three 
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areas.  First,  there  are  photoelastic  methods,  such  as  those  used  by  Cornell  [3]  to  augment  his 
theoretical  work  mentioned  earlier,  or  Hahn  [15],  Kutscha  [16]  and  others.  Hahn's  work  presents 
two  interesting  conclusions,  which  have  not  been  investigated  in  the  literature.  First,  he  notes 
inelastic  behavior  in  the  adherend,  which,  as  previously  discussed,  has  not  been  a  focus  of 
analytical  efforts.  Second,  Hahn  notes  photoelastic  behavior  in  the  bond  layer,  particularly  at 
the  end  of  lap,  which  he  interprets  as  a  variation  in  stress  across  the  adhesive  layer.  This  is 
contrary  to  the  assumption  of  infinitessimal  thickness  and  continuous  tractions  from  adherend 
to  adherend  as  used  in  classical  solution.  The  existing  numerical  techniques  and  analytical 
methods  fail  to  show  this  variation  across  the  bond  line.  This,  then,  is  a  major  objective  of  the 
research  reported  herein. 

The  second  area  is  the  modeling  of  the  real  metal-adhesive-metal  joint  by  rubbers  of 
moduli  appropriate  to  the  modeling  laws  given  above.  Early  work  was  done  by  NASA  [17]  and 
more  recent  work  has  been  done  by  Adams,  et  a  I .  [13,18].  The  primary  problem  with  both  the 
photoelastic  methods  and  rubber  models  is  that  only  a  qualitative  idea  of  the  stress  distribution 
can  be  determined.  Experimental  models  in  general,  and  the  rubber  model  in  particular,  have 
three  additional  drawbacks.  First  it  has  not  been  shown  that  the  behavior  of  a  one  half  inch 
or  greater  model  adhesive  layer  is  a  good  approximation  to  the  real  problem.  Second,  although 
the  Young's  modulus  ratio  or  the  shear  modulus  ratio  are  similar  for  the  rubber  model  and  the 
real  problem,  the  Poisson  ratios  are  greatly  different.  It  has  not  been  shown  that  this 
difference  has  a  negligible  effect  on  the  results.  Third,  the  rubber  models  are  actually  a  five 
layer  system  with  stiff-rubber/bond/foam-rubber/bond/stiff-rubber.  Although  the  bond  has 
similar  properties,  it  is  not  clear  that  the  five  layer  system  closely  approximates  the  real 
problem. 
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The  final  area  of  experimental  work  is  that  of  destructive  testing  of  actual  adhesive 
bonded  joints  to  determine  allowable  loading,  mechanisms  of  fracture  and  other  data  required 
for  safe  design  of  joints.  The  Forest  Products  Laboratory  Report  [2]  and  research  reported  by 
Wang,  et  a  I .  [4],  Dukes  and  Bryant  [19],  and  Crimes  [20],  are  typical  examples  of  destructive 
testing.  Although  significant  observations  can  be  made  by  direct  experimentation,  and  the 
ultimate  verification  of  any  analytical  method  is  by  direct  experimentation,  this  technique  does 
not  lead  to  understanding  of  the  general  case.  That  is  to  say,  while  Dukes  and  Bryant  have 
established  a  relationship  between  adhesive  thickness  and  reduced  strength,  it  is  known  to  be 
valid  only  in  conjunction  with  the  geometry  and  materials  variables  used.  If  an  analytical  or 
numerical  method  can  be  perfected  which  can  predict  failure  from  the  computed  stress  field, 
and  if  it  can  be  verified  by  destructive  testing,  then,  and  only  then,  will  there  be  a  high  degree 
of  confidence  in  the  design  of  adhesive  joints,  without  a  costly  testing  program  of  actual 
geometries  and  materials.  The  research  reported  herein  is  hoped  to  be  a  step  toward  such  a 
method. 

Section  III:  Finite  Element  Method 

The  FEM  is  an  extremely  versatile  and  powerful  application  of  the  variational  principles 
of  mechanics  by  which  a  complex  shaped  continuum,  such  as  the  single  lap  shear  joint,  may  be 
divided  into  a  finite  number  of  discrete  elements.  After  dividing  the  continuum,  displacement 
or  stress  fields,  or  both,  are  assumed  within  each  element.  The  equations  that  result  from  the 
application  of  the  appropriate  variational  principle  are  simultaneous  linear  algebraic  equations 
which  have  generalized  displacements,  as  unknowns  to  be  determined  [21]. 
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In  the  research  reported  herein,  three  different  elements  were  used  (not  counting  the 
super  element  used  at  the  crack  tip  of  the  cases  that  contained  a  crack.  See  Section  IV.  Two 
elements,  an  isoparametric  triangular  element,  and  an  isoparametric  quadrilateral  element,  are 
derived  from  the  principle  of  minimum  potential  energy  and  are  briefly  outlined  in  section 
Ill.i.  below.  The  remaining  element  is  a  rectangular  hybrid  element  [22],  based  on  the 
complimentary  energy  principle  outlined  in  Section  Ill.ii.  below.  The  reason  for  using  both  the 
isoparametric  quadrilateral  element  and  the  assumed  stress  hybrid  rectangular  element  in  the 
same  problem  is  primarily  one  of  accuracy.  The  rectangular  hybrid  element  is  as  accurate  as 
an  eight  node  isoparametric  quadrilateral  element,  but,  having  fewer  degrees  of  freedom,  is  far 
more  efficient  with  respect  to  computation  time.  The  assumed  stress  hybrid  method  does  not 
require  a  rectangular  shape,  however,  if  a  rectangular  shape  is  used,  integrations  along  the 
boundary  can  be  done  before  programming,  which  greatly  increases  efficiency.  A  general 
quadrilateral  shape  using  the  hybrid  model  requires  numerical  integration  along  each  boundary 
segment  which  increases  computation  time.  This  problem  was  avoided  in  the  research  reported 
herein  by  arranging  the  mesh  such  that  all  general  quadrilater  shapes  are  only  in  less  critical 
areas.  Therefore,  the  more  efficient  (one  nine  point  numerical  integral  vice  eight  five  point 
numerical  line  integrals),  but  less  accurate,  isoparametric  element  was  used  in  these  areas,  rather 
than  a  general  quadrilateral  assumed  stress  hybrid  element. 
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This  means  that  the  entire  problem  reduces  to  the  form 

Kq  =  Q 

where 

K  =  the  assembled  global  stiffness  of  the  entire  structure 

q  =  generalized  nodal  displacements 

Q  =  equivalent  nodal  forces 

All  three  elements  use  a  linear  displacement  interpolation  function  which  assures 
compatibility  along  interelement  boundaries.  Therefore,  it  is  appropriate  and,  as  it  will  be 
shown,  advantageous  to  use  elements  formulated  from  different  variational  principles. 

i.  Isoparametric  Elements 

The  isoparametric  elements  are  based  on  the  assumption  of  a  displacement  field 
continuous  over  the  entire  solid.  For  this  reason,  they  are  often  referred  to  as  compatible 
models. 

The  formulation  of  the  isoparametric  triangular  element  or  the  isoparametric 
quadrilateral  element  can  be  found  in  any  introductory  texts  on  the  FEM  [8,91  The 
isoparametric  triangular  element  is  first  transformed  to  a  parent  element,  as  shown  in  Figure  4, 
in  triangular  or  area  coordinates  such  that 

3  3 

x  =  2  pi  xi  y =  2  pi  y* 
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where 


P|-tl 


and 

£1  ♦  £2  ♦  $3  =  ' 

Similarly  for  the  quadrilateral  element,  there  is  a  transformation  to  a  parent  element,  as 
shown  in  Figure  5,  such  that 

*  =  I  Pi  x,  y  -  2  Pi  y, 


where 


and 


Pi  -   1/4  (I  ♦  &)  (1  ♦  W.) 


*■  .    2p.  =  l 
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The  principle  of  minimum  potential  energy  for  a  plane  elasticity  problem  may  be  stated 
as  the  vanishing  of  the  variation  of  the  total  potential  energy  functional,  7Tp,  which  can  be 
written  for  each  element  as 

ftp  - 1  /A  (0.5  €T  E   e  -   u  F)  dA   -   fsa.     u  T  ds 


where 


€  =  strain  vector 


E  -  elastic  constants 


A  ■  area  of  the  element 


F  ■  prescribed  body  forces 

T  »  prescribed  surface  traction 

U  ■  generalized  displacements 

s  »  boundary  curve  of  element 

$a  »  portion  of  s  over  which  surface  tractions  are  prescribed 

t  ■  thickness  (assumed  unity) 
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But  U  can  be  written  as: 


u  =  Rq 


where 


q  =  unknown  nodal  displacements, 


R  =      |Rj  R2 Rn(    n=3  for  triangular  element,  n=4  for  quadrilateral  element 


and 


R.= 


Pi      o 
0         P. 


therefore 


e  =  Du 


where 


D  is  the  differential  operator 


(blbx) 

0 

(blby) 


0 
O/dy) 

(blbx) 


> 


e  can  be  written 


e  =  Bq 


where 


?B,B2 Bn 


Bi 


(bVJbx)        0 

0  (b?Jby)    r 

(bPJby)        (bVJbx) 


defining  the  element  stiffness  matrix 

k  =  /A  b'  E  B  dA 


and  the  equivalent  nodal  force 

Q  =  /A    RTFdA  +  /S(_RTds 


The  interpolation  functions,  however,  are  not  in  terms  of  x  and  y  so  a  modification  must 


be  made  to  the  integration  to  account  for  that  fact.  That  is 

T  -IT 


k  =  /AB  [J  Y  E  B[J]  dA 
Q  =  /A   RFdA  +  Z,,  RTds 


where 


dA 


dA 


&  (x,y) 
&(x,y) 


d£i  d%2       f°r  tne  triangular  element 


2Xfo> 


d£di]    for  the  quadrilateral  element 


The  integration  of  the  triangular  element  simplifies  significantly  due  to  the  property  of 
integration  in  triangular  coordinates 

where 

&(x,y) 
=  2A 

and 


A  =  Area  of  element 


then 


m  !  n  !  s ! 
2    3  (rmn+s+2)! 


The  integration  of  the  terms  of  k  for  the  quadrilateral  element  are    done  numerically 
using  a  nine  point  Gaussian  quadrature  of  the  form  (bi-linear  three  point  quadrature) 


/B  Mx.  y)  dA  -  f!    W,W,  f(x,.y) 


where 


Wj.Wj  =•  weighting  factors 

f\x,,y)  ■  value  of  the  function  /"(x.y)  evaluated  at  {x,y) 
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ii.  Rectangular  Hybrid  Element 

The  assumed  stress  hybrid  model  is  based  on  the  complementary  energy  principle.   That 
Drinciple  can  be  written  by  the  following  functional  for  any  element  [22]: 

7TC  -  t  /R  0.5  <r  C  <rdA  -    t/AuF   dA  -     /8u  U  T    ds 

where 

<r  =  stress  component  matrix 

C  =  elastic  compliance  matrix 

U  ■  prescribed  displacements  over  sy 

u  "  displacement  over  s 

T  =  boundary  traction  over  s 

A  =  area  of  element 

F  =  body  forces 

Sy  =  boundary  of  A  over  which  displacements  are  prescribed 

s  -  boundary  of  A 

t  =  thickness  (assumed  unity) 
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The  variational  principle  then  states  that  of  all  the  admissable  stresses  and  displacements, 
;he  actual  solution,  that  is  the  one  that  satisfies  compatibility,  equilibrium  and  boundary 
conditions  in  A,  is  distinguished  by  the  stationary  value  of  7TC  with  respect  to  variations  of  <r  in 
A  and  u,  over  s  and  su. 

The  extremum  of  the  functional  is  also  subject  to  the  constraint  conditions 

V<r+    F    =    0 
<rv  =   T 

where 

v  ■   unit  normal  vector  to  s 

Simply  stated,  the  admissible  stresses  must  be  continuous  in  A  and  satisfy  the  equilibrium 
equations,  and  be  compatible  with  prescribed  boundary  tractions.  Admissible  displacements 
must  satisfy  the  prescribed  boundary  displacements  over  sy. 

In  the  hybrid  model,  the  stress  in  the  element  and  the  displacement  along  the  boundary 
are  assumed  separately.  A  typical  element  is  shown  in  Figure  6.  For  such  an  element,  the 
stresses  are  given  as 

<r  =  P  0 

where 

P  ■  unknown  stress  parameters  satisfying  the  equilibrium  equations. 
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The  only  requirement  on  the  number  of  /3s  is  that  there  must  be  a  sufficient  number  to 
avoid  unstable  kinematic  modes.   Let 

n  ■  number  of  nodal  degrees  of  freedom 

I  -  number  of  rigid  body  modes 

» 

m  ■  number  of  /3s 

A  necessary  condition  for  a  stable  solution  is  that  m  >  n-1,  which,  for  the  case  of  the 
rectangular  element  with  four  nodes  with  two  degrees  of  each,  means  there  must  be  five  or 
more  /Ss.  This  is,  however,  only  a  necessary  condition  and  no  sufficient  condition  is  known  to 
exist.  The  research  reported  herein  uses  a  five  /S  system  the  accuracy  of  which  is  indicated  by 
Plan  [211 


Then 


r 


o      o 


0         10         0 


The  boundary  tractions  can  be  expressed  as: 


T  =  R/3 


.                   v. 

"12 

T 

*I2 

"41 

Tu 

*41 

\            • 

^ 


x  V 


0         0         10         0 


but  for  a  rectangular  element  T  reduces  to 


0-10 


-1 


0         -x 


1  0         0  y         0 

0  0  10         0 


0  1 


0 
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0  1  0         0  x 

-1  0  0         -y         0 

0         0-100 


In  this  hybrid  model,  the  interelement  boundary  displacement  function  Uj  is  interpolated 
in  terms  of  the  nodal  displacements  q.  That  is 

u  =  L  q 


where 


(1-x/a) 

0 

x/a 

0 

0 

0 

0 

0 

0 

(1-x/a) 

0 

x/a 

0 

0 

0 

0 

0 

0 

(1-y/b) 

0 

y/b 

0 

0 

0 

0 

0 

0 

(1-y/b) 

0 

y/b 

0 

0 

0 

0 

0 

0 

x/a 

0 

(1-x/a) 

0 

0 

0 

0 

0 

0 

x/a 

0 

(1-x/a) 

(1-y/b) 

0 

0 

0 

0 

0 

y/b 

0 

0 

(1-y/b) 

0 

0 

0 

0 

0 

y/b 

Defining"  the  following  matrices 


H  =  /A  PT  C  P  dA 


G  =  /^  RT 


dA  R  L  ds 
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H  is  determined  by  a  nine  point  Gaussian  quadrature  over  the  area.   G  is  determined  by 
five  point  quadrature  along  each  applicable  boundary  segment.   The  variational  functional 
fc  can  then  be  written 

TCC  =  1/2  0T  H/3  -   0TGq 

The  stationary  condition  of  the  functional  with  respect  to  0  and  q  gives 

H/3  -  Gq  =  0 

Solving  for  /S  in  terms  of  q  gives 

0  =  H-'Gq 

If  the  following  is  observed 

k  =  GH-'G 

Then  the  variational  functional  can  be  written  in  a  form  identical  to  that  of  the 
soparametric  elements  in  the  previous  section,  i.e., 

7TC  -  1/2  q'kq  -q  Q 

The  individual  element  stiffness  matrices  can  then  be  assembled  systematically  node  by 
lode  into  a  global  stiffness  matrix.  This  forms  a  set  of  matrix  equations,  which,  when 
:ombined  with  the  boundary  conditions  of  the  problem,  yield  the  solution  to  the  problem.  Both 
he  isoparametric  elements  and  the  assumed  stress  hybrid  model  use  linear  interpolation  of 
nodel  displacements,  so  that  the  resulting  solution  to  the  set  of  assembled  matrix  equations  will 
iave  compatible  displacements  over  the  whole  domain. 
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ection  IV:  Hybrid  Stress  Finite  Element  Analysis  in  Fracture  Mechanics  of  Adhesive  Joints 

The  study  of  fracture  of  adhesive  bonded  joints  has  become  of  great  commercial  interest 
s  the  use  of  such  joints  becomes  more  widespread.  The  study  of  cracks  in  adhesive  joints 
ippears  broken  into  two  general  areas.  First,  there  are  the  primarily  experimental  studies,  such 
is  those  of  Mostovoy  and  Ripling  [24]  or  Jermian  and  Ventrice  [25].  Second,  there  is  the  study 
bf  fracture  by  use  of  the  FEM  by  both  conventional  displacement  models,  such  as  Wilson  [26], 
Trantina  [26,27],  Anderson,  et  al  .  [28],  and  Williams  [29],  as  well  as  hybrid  stress  models 
which  include  a  singularity  in  the  assumed  stress  field  [23,30-32]. 

Although  the  experimental  determination  of  the  strain  energy  release  rate  in  a  cracked 
monolithic  material  is  widely  used,  there  is  some  question  as  to  its  applicability,  particularly  to 
the  lap  shear  configuration.  For  the  symmetric  case,  such  as  the  double  edge  notch  specimen, 
the  results  for  Gj  [24]  compare  very  closely  with  those  calculated  by  an  exact  conformal 
mapping  technique  by  Bowie  [33].  Gj  and  kj  are  related  for  a  monolithic  material  through  the 
expression: 

V  (k  *  1)  kj 

S,  -  --" 

8G 

where 

G  ■  shear  modulus 

K  »  3-4y  for  plane  strain  and  (3-y)/(l-«r)  for  plane  stress 
v  ■  Poisson's  ratio 
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Current,  yet  to  be  published  work  by  McGarry,  Wang  and  Mandeli  at  M.I.T.  shows  that 
the  above  expression  also  applies  to  adhesive  joints  if  the  shear  modulus  and  Poisson  ratio  of 
the  adhesive  are  used.  This  is  true  because  the  characteristic  length  of  the  singularity  at  the 
crack  tip  is  much,  much  less  than  the  characteristic  length  of  the  joint  which  was  found  to  be 
the  adhesive  thickness.  The  basis  of  this  experimental  approach  is  to  determine  dC/da,  or  the 
change  in  compliance  with  respect  to  crack  length.  In  the  double  cantilever  beam  case,  simple 
beam  theory  gives  a  very  good  approximation  of  the  compliance  as  the  crack  length  changes, 
so  it  is  not  unexpected  that  the  results  are  reasonable  [24].  In  the  lap  shear  configuration,  the 
determination  of  the  compliance  is  a  non-trivial  problem,  due  to  the  complex  geometry,  its 
unknown  relationship  with  rotation,  and  differential  straining  of  the  joint.  The  strain  energy 
release  rate  approach  was  used  by  Jermian  and  Ventrice  [25],  but  how  they  arrived  at  their 
expression  for  dC/da  is  not  clear  from  their  reported  research. 

The  normal  displacement  model  FEM  has  a  great  deal  of  difficulty  handling  the  crack 
tip  region  in  a  fracture  problem,  even  in  the  cracked  monolithic  material.  Even  with  an 
extremely  fine  mesh  in  the  vicinity  of  the  crack  tip,  such  as  that  used  by  Trantina  [27],  the 
estimate  of  Gj  is  only  257.  of  that  computed  by  Bowie  for  the  symmetric  problem.  Anderson,  ct 
al.  [28]  and  Williams  [29]  have  applied  relatively  large  constant  strain  elements  to  the 
interfacial  fracture  problem.  Due  to  the  type  of  element  (constant  strain),  and  the  relative 
coarseness  of  the  mesh  (.002  inch)  for  the  type  of  problem,  it  is  not  surprising  that  their 
solution  shows  sharp  oscillations  about  the  analytic  solution  [34]. 

The  first  work  in  the  FEM  to  incorporate  stress  singularities  in  the  assumed  stress  hybrid 
model  contained  only  a  term  in  1/-/F  [23,30].    The  later  work  includes  higher  order  singular 
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terms  to  more  accurately  account  for  the  stress  singularity  at  the  tip  of  the  crack  [31,323.  It  is 
this  latter  approach  which  is  used  as  the  basis  for  the  super  element  used  at  the  crack  tip  in  the 
research  reported  herein. 

The  formulation  of  the  super  element  uses  the  complex  variable  stress  function 
formulation  of  Muskhelishvili  [37].  The  H  and  G  matrix  are  then  formed  by  line  integrations 
in  the  complex  plane  along  each  segment,  shown  in  Figure  7.  Although  the  symmetic  five  node 
element  has  been  used  and  gives  good  agreement  with  the  exact  solution,  the  complete  nine 
node  super  element  has  never  been  perfected  prior  to  the  research  reported  herein.  During  the 
course  of  this  study,  it  was  found  that  the  last  two  rows  of  the  stiffness  matrix  published  by 
Lasry  [31]  are  in  error  and  some  of  the  formulations  in  the  paper  by  Tony  [32]  are  also 
incorrect.  The  correct  stiffness  matrix  for  the  nine  node  super  element  is  included  as  Table  3. 
The  stiffness  matrix  given  by  Lasry  yields  non-symmetric  displacements  for  the  problem  shown 
in  Figure  8.  It  should  be  noted  here  that  the  Aeroelastics  Laboratory  at  MIT  has  corrected 
their  solution  of  the  stiffness  matrix. 

The  stress  intensity  factors  are  shown  [32]  to  be: 

k,    -   720, 

The  research  reported  herein  uses  N  =  9,  which  meets  the  previously  stated  necessary 
condition  to  avoid  unstable  kinematic  modes.  The  assumed  stress  hybrid  super  element  is  also 
a  stiffness  method,  and  the  resulting  element  stiffness  matrix  can  be  assembled  as  with  the 
element  stiffness  matrices  of  the  other  elements  to  form  the  matrix  equations.    The  super 
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element  also  uses  a  linear  interpolation  function  between  each  two  nodes  which  ensures  inter- 
element  compatibility  of  the  displacement  solution  with  the  other  elements  discussed  above. 

Section  V:  Results  of  Analysis  of  Joint  Without  Crack 

The  single  most  important  factor  to  examine  when  solving  a  problem  by  a  numerical 
technique  is  the  stability  of  the  solution  with  respect  to  changes  in  the  characteristic  parameters 
of  the  method.  That  is  to  say,  although  convergence  proofs  do  exist  for  the  FEM  [35],  it  can 
only  be  said  that  the  finite  element  solution  approaches  the  exact  solution  as  the  mesh  size 
becomes  "small."  The  first  task  in  a  numerical  solution  scheme  is  to  determine  what  "small"  is 
for  the  specific  problem  at  hand.  In  a  problem  for  which  an  exact  analytic  solution  exists,  this 
determination  can  be  made  by  selecting  the  mesh  arrangement  that  is  acceptably  close  to  the 
analytic  solution.  In  the  adhesive  bonded  lap  shear  problem,  no  such  exact  solution  exists. 
The  fact  that  no  exact  solution  does  exist  requires  that  a  stable  solution  be  found. 

Determination  of  a  stable  mesh  size  must  be  done  systematically  for  two  reasons.  First, 
the  computation  time  required  to  solve  the  problem  is  proportional  to  the  square  of  the  band 
width  of  the  matrix  equations  and  linearly  proportional  to  the  total  number  of  degrees  of 
freedom.  That  is 

t  =»c(  n2  +  d.o.f. ) 

where 

n  =  band  width  of  matrix  equations 
d.o.f.  *  total  degrees  of  freedom 
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Therefore,  the  decrease  of  mesh  dimensions  has  a  highly  non-linear  effect  on 
computation  time.  Second,  a  mesh  that  is  too  small  may  give  an  inaccurate  solution  to  the 
problem  due  to  round  off  error  in  computation.  Roundoff  error  is  greatly  reduced  in  the 
research  reported  herein  by  doing  all  computation  in  a  double  precision  mode. 

Figure  9  shows  a  plot  of  the  peel  stress,<ry,  at  the  end  of  lap  along  the  interface  between 
adherend  and  adhesive  for  various  meshes.  The  stability  was  checked  in  locations  other  than 
along  the  interface  and  for  stresses  other  than  the  peel  stress,  but  the  results  with  respect  to 
solution  stability  were  the  same.  The  latter  is  not  surprising  when  one  considers  that  it  is  the 
displacement  field  which  is  converging  on  the  correct  solution.  Because  all  stresses  are  a 
function  of  nodal  displacement  only,  all  stresses  should  converge  in  a  similar  manner. 

Figure  10  shows  the  mesh  arrangement  for  the  critical  region  in  the  vicinity  of  the  end  of 
lap.  Of  primary  interest  is  the  use  of  eight  elements  across  the  adhesive  layer.  The  large 
number  of  elements  across  the  adhesive  thickness  provides  for  the  definition  of  the  variation 
of  stress  across  the  thickness.  This  has  not  been  reported  previously,  either  analytically  or 
numerically,  although  it  has  been  stated  qualitatively  by  Hahn  [15]  from  his  experimental 
research.  Figure  12  shows  the  variation  in  normal  or  peel  stress  across  the  adhesive  thickness  at 
various  distances  along  the  lap.  Figure  13  is  a  similar  plot  for  shear  stress. 

The  strain  across  the  interface  physically  must  be  continuous,  since  both  materials 
(adherend  and  adhesive)  are  continuously  connected.  Because  the  strain  is  continuous, 
longitudinal  stress  will  be  different  in  the  two  materials  due  to  different  elastic  properties.  To 
satisfy  compatibility,  however,  peel  stress  and  shear  stress  (rxy)  must  be  continuous.  It  is 
significant  to  note  that  although  it  is  possible  to  introduce  a  constraint  condition   into  the 
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formulation  of  the  problem,  the  imposition  of  strict  stress  continuity  in  cry  and  rxy  across  the 
adherend  adhesive  interface  was  not  done.  Therefore,  in  determining  the  value  of  <ry  or.rxy  at 
the  interface,  an  extrapolation  must  be  done,  both  from  the  adhesive  side  and  the  adherend 
side.  If  the  two  values  agree,  there  is  a  high  degree  of  confidence  in  the  actual  value  along  the 
interface.  Such  extrapolations  from  both  sides  of  the  interface  were  done  for  all  values 
presented  on  the  interface,  and  it  was  found  that  the  two  extrapolations  agreed  within  plotting 
accuracy.   Thus,  there  is  a  strong  case  for  the  reliability  of  stress  values  along  the  interface. 

Figure  12  and  13  present  the  variation  in  stress  across  the  adhesive  thickness  for  only  one 
modulus  ratio  of  the  nine  ratios  used.  It  was  found  that  the  same  trend  appears  throughout: 
the  gradient  in  stress  is  greater  near  the  end  of  the  lap.  The  lap  is  highly  non-linear  in  stress 
near  the  upper  adherend  and  becomes  roughly  linear  within  about  one  adhesive  thickness  and 
nearly  constant  across  the  thickness  in  about  five  adhesive  thicknesses.  The  non-linear 
behavior  near  the  end  of  the  lap  is  much  more  pronounced  for  low  modulus  ratios. 

As  discussed  in  Section  II,  the  analytical  and  numerical  literature  on  adhesive  bonded  lap 
joints  does  not  address  longitudinal  stress  in  the  adherend.  This  is  a  serious  shortcoming  in 
light  of  the  experimental  evidence  [4,15]  of  failures  occuring  in  the  adherend  due  to  local 
bending  at  the  end  of  the  lap.  Figure  14  shows  the  longitudinal  stress  in  the  adherend  at  the 
interface;  since  the  strain  is  continuous,  the  longitudinal  stress  in  the  adhesive  is  much  lower., 
As  shown  in  Figure  14,  however,  the  stress  concentration  may  be  sufficiently  high  to  cause 
yielding,  or  failure,  if  a  low  strength  adherend  such  as  1100  aluminum  is  used. 

The  longitudinal  stress  concentration  reduces  somewhat  as  the  modulus  ratio  increases. 
This  reduction  is  to  be  expected  as  differential  straining  becomes  more  important  relative  to 


rotation  and  as  the  adhesive  becomes  more  flexible.  The  reduction  in  longitudinal  stress  with 
change  in  moduli  ratio  is  small,  however,  varying  only  from  4.8  for  a  modulus  ratio  of  one  to 
about  4.1  for  a  modulus  ratio  of  one  thousand. 

For  a  fixed  modulus  ratio,  the  longitudinal  stress  concentration  increases  with  increasing 
adhesive  thickness.  This  is  due  to  the  increased  rotation  caused  by  increased  eccentricity  of  the 
applied  loading.  The  magnitude  of  the  longitudinal  stress  concentration  varies  from  4.2  for  an 
adhesive  thickness  of  2.5  mils  to  4.9  for  an  adhesive  thickness  of  20  mils  (fixed  modulus  ratio 
of  20). 

Figures  15  through  20  give  additional  data  on  the  variation  in  stress  from  one  interface 
to  the  other.  Previous  solutions  to  the  problem  have  indicated  that  the  stresses  are  symmetric 
about  the  center  of  the  lap.  Figures  15  and  16  show  the  error  in  that  assumption.  In  both  the 
moduli  ratios  of  commercial  importance  in  adhesive  joints  (20  and  200),  shown  in  Figures  15 
and  16,  the  fillet  shape  at  the  end  of  the  lap  has  greatly  relaxed  the  stress  in  the  vicinity  of 
point  A  (Figure  15)  compared  to  the  stress  in  the  vicinity  of  point  B  (Figure  16).  Modulus  ratio 
2  shows  a  behavior  more  nearly  like  a  homogeneous  material,  i.e.,  it  is  very  strongly  influenced 
by  geometric  discontinuities.  That  is,  modulus  ratio  two  shows  an  inflexion  indicating  a 
relaxation  of  stresses  similar  to  modulus  ratios  20  and  200,  but  upon  continuing  toward  the 
geometric  discontinuity,  the  singular  behavior  at  the  discontinuity  overwhelms  the  relaxation 
trend. 

Figure  17  and  Figure  18  are  the  same  data  as  Figures  15  and  16,  but  "normalized"  by  the 
value  of  ay  at  point  B^.  Point  B^  is  the  value  of  stress  at  the  centroid  of  the  closest  element  to 
point  B,  extrapolated  to  the  interface  by  the  method  discussed  above.   The  distance  from  point 
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B,  which  is  the  actual  end  of  lap,  to  point  B/  is  0.1  adhesive  thicknesses.  Although  the  general 
shape  of  the  curves  is  not  affected,  Figures  17  and  18  do  make  it  easier  to  see  the  distance  to 
which  the  perturbation  is  observed  in  the  joint.  Much  higher  absolute  values  of  peel  stress  are 
found  for  low  moduli  ratios,  as  shown  in  Figures  15  and  16,  but  Figures  17  and  18  show  how 
rapidly  these  decay.  Perturbation  distance  is  of  the  order  of  only  one  adhesive  thickness  for 
modulus  ratio  2,  and  increases  for  the  more  flexible  modulus  ratios  to  about  5  thicknesses  for 
ratio  200. 

Figures  19  and  20  are  the  shear  stresses  along  the  interface  normalized  by  the  value  at 
point  B^.  Here,  again,  the  fillet  shape  at  point  A  has  virtually  eliminated  the  stress 
concentration  in  the  higher  moduli  ratios.  The  perturbation  distance  for  shear  is  greater  than 
that  for  peel,  as  would  be  expected  from  the  eccentric  loading  of  the  single  lap  joint. 

Figures  21  and  22  give  the  peeling  stress  and  shear  stress  at  the  midplane  of  the  bond. 
As  seen  in  Figures  2  and  3,  the  midplane  solution  is  very  close  to  the  Goland  and  Reissner 
solution  in  the  critical  end  of  the  lap  region.  As  the  earlier  figures  have  shown,  the  maximum 
stress  occurs  at  the  adherend-adhesive  interface,  and  the  midplane  is  a  conservative  location  of 
stress  concentration.  The  numerical  technique  of  Wooley  and  Carver  [10]  loses  the  peak  stress 
of  point  B'  due  to  the  use  of  only  two  constant  strain  elements  across  the  thickness.  If  the 
values  for  stress  obtained  in  the  research  herein  are  arithmetically  averaged  through  the 
thickness,  they  are  very  close  to  the  other  numerical  solutions  (which  internally  average  due  to 
the  nature  of  the  constant  strain  elements),  such  as  those  of  Wooley  and  Carver,  which  lends 
weight  to  the  reliability  of  the  solution  reported  herein. 

Figure  23  is  included  to  show  the  variation  in  the  maximum  peel  stress  and  maximum 
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ihear  stress  at  point  B^  as  the  modulus  ratio  changes  over  three  orders  of  magnitude.  The 
-naximum  stress  decreases  monotonically  as  modulus  ratio  increases.  This  is  because  the  more 
flexible  adhesive  layers  can  distribute  the  effect  of  geometric  discontinuities  over  a  greater 
distance  as  shown  by  larger  perturbation  distances  for  the  higher  moduli  ratios. 

Figure  24  shows  the  variation  in  peel  stress  along  the  interface  in  the  vicinity  of  point  B 
for  various  adhesive  thicknesses.  From  the  Figure,  it  is  seen  that  the  maximum  change  in 
slope  for  all  thicknesses  occurs  at  roughly  the  same  location,  which  supports  the  earlier  finding, 
i.e.,  that  perturbation  distance  is  primarily  a  function  of  the  moduli  ratio.  The  results  of 
Figure  24  must  be  interpreted  carefully  for  the  following  reasons.  First,  the  analysis  scheme 
used  does  not  take  into  consideration  thermal  elastic  effects  or  residual  stress,  both  of  which 
become  more  important  for  thicker  adhesive  layers.  That  means  that  the  fact  that  the  loading 
condition  does  not  produce  a  large  stress  concentration  in  the  joint  may  be  unimportant,  if  the 
joint  is  at  the  point  of  incipient  failure  due  to  shrinkage-induced  residual  stresses.  Second, 
failure  is  statistical  in  nature,  and  the  larger  the  joint,  the  more  likely  a  failure  producing  flaw 
is  to  exist  [19].  Third,  a  joint  with  a  thicker  adhesive  layer  will  have  less  dimensional  stability 
due  to  the  lower  modulus  and,  therefore,  larger  strains  in  the  adhesive.  Finally,  as  shown 
earlier,  the  longitudinal  stress  in  the  adherend  is  increased  for  greater  adhesive  thicknesses,  and 
may  induce  plastic  deformation,  or  even  failure,  in  the  adherend.  The  seriousness  of  these 
thickness  effects  can  be  appreciated  when  one  considers  that  the  range  of  adhesive  thicknesses 
normally  used  fall  within  the  range  of  three  to  ten  mils. 
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Section  VI:  Results  of  Analysis  of  Joint  with  Short  Edge  Crack 

As  discussed  in  Section  V,  the  first  step  in  the  numerical  solution  of  a  problem  for  which 
no  analytic  solution  exists  is  to  determine  stability  of  the  solution.  The  determination  of  a 
stable  solution  when  using-  an  assumed  stress  hybrid  super  element  is  more  complicated  due  to 
the  nature  of  the  element.  Unlike  the  conventional  displacement  model  of  the  FEM,  where 
smaller  elements  yield  a  closer  approximation  to  the  exact  solution  (not  countng  round  off 
error),  the  super  element,  being  a  numerically  exact  solution  including  singular  and  higher 
order  terms  [33],  yields  a  better  solution  for  a  larger  element.  Because  the  crack  tip  super 
element  is  surrounded  by  normal  elements,  however,  the  overall  quality  of  the  stress  field 
determination  is  a  trade  off  between  a  super  element  that  is  as  large  as  possible,  and 
surrounding  elements  that  are  as  small  as  practical. 

The  super  element  must  be  sufficiently  large  to  contain  the  singular  behavior.  Because 
the  stress  gradient  is  still  very  large  at  the  interelement  boundary  between  the  super  element 
and  the  surrounding  elements,  a  study  of  the  mesh  size  in  the  vicinity  of  the  crack  tip  is 
necessary  to  determine  a  stable  solution.  The  research  reported  herein  did  not  systematically 
determine  the  optimum  mesh  arrangement.  As  seen  in  Figure  25,  however,  the  I/V7  singular 
behavior  of  stress  does  end  within  the  superelement.  In  addition,  a  vast  amount  of  work  has 
been  done  on  the  fracture  problem  of  adhesive  joints  by  Dr.  Wang,  M.I.T.  Materials  Science 
and  Engineering  Department,  using  the  five  node  symmetric , super  element.  In  the  course  of 
his  work,  optimum  mesh  arrangements  were  determined  for  a  similar  adhesive  joint  system, 
and,  therefore,  the  mesh  arrangement  shown  in  Figure  11,  which  was  used  in  the  work  reported 
here,  was  not  a  purely  arbitrary  choice.    The  mesh  chosen  conformed  to  the  optimum  mesh 
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arrangement  determined  by  Dr.  Wang's  work,  and,  therefore,  the  solution  reported  herein  is 
estimated  to  be  within  one  half  percent  of  the  converged  solution.  Because  no  solution  exists  to 
the  problem  of  an  adhesive  bonded  lap  joint  with  an  edge  crack,  the  results  reported  herein  are 
a  significant  initial  input  to  the  continuing  research  on  this  problem. 

One  interesting  result  arising  from  the  analysis  of  the  adhesive  bonded  lap  joint  with  a 
short  edge  crack  is  shown  in  Figure  26.  This  figure  is  a  plot  of  log  stress  intensity  factor  kj 
(opening  mode)  and  kll  (inplane  shear  mode)  vs  log  modulus  ratio,  where  these  ks  are  related 
to  the  conventional  nomenclatures  used  by  NASA  and  ASTM  by  the  following: 

K,   -   Vi   k, 

Kn-V?  kH 

Figure  26  is  significant  in  that  at  least  for  moduli  ratios  between  2  and  200,  both  the 
opening  mode  stress  intensity  factor  and  the  inplane  shear  mode  stress  intensity  factor  plot  as 
straight  lines  in  log-log.  The  slopes  of  these  lines  are  0.49  for  kj  and  0.46  for  kn.  This 
implies  that 

k,   -  C,   (E/E/<9 
k„-C2  (E/E,)M6 

where  C,  and  C2  are  constants  dependent  on  the  nominal  load,  crack  length  and  other 
geometric  variables  of  the  problem.  Such  a  relation  must  be  systematically  investigated  further 
to  determine  if  the  power  law  relating  the  modulus  ratio  to  the  stress  intensity  factor  is  itself  a 
function  of  the  geometric  variables. 

A  comparison  of  the  solution  to  the  lap  shear  problem  with  and  without  a  short  edge 
crack  permits  the  following  observations: 


1.  Interfacial  stresses  at  the  end  of  the  lap  without  a  crack  are  identical  to  those 
found  in  Section  V  for  the  uncracked  solution. 

2.  Longitudinal  stress  concentrations  in  the  adherend  along  the  interface  are  the 
same  as  the  uncracked  case  except  that  the  location  of  the  maximum 
longitudinal  stress  moves  toward  the  center  of  lap  by  approximately  one  half 
the  crack  length. 

3.  Far  field  stresses  in  the  adherends  are  identical  for  both  cracked  and  uncracked 
solutions. 

4.  Stresses  in  the  crack  flank  are  relaxed,  from  those  found  in  the  end  of  lap  of 
the  uncracked  case.  (Figures  15-20  vs.  Figures  27-29) 

The  first  three  points  indicate  that  the  presence  of  a  short  crack  does  not  affect  the 
stiffness  of  the  joint  as  a  whole.  If  the  stiffness  were  significantly  affected,  there  would  be  an 
increased  longitudinal  stress  from  the  additional  bending  required  to  accommodate  additional 
rotation,  or  the  far  field  stresses  would  begin  to  change,  or  both. 

By  comparing  Figures  27  and  2S  with  Figures  16  and  20,  respectively,  the  drastic  effect 
on  the  interfacial  stresses  along  the  upper  interface  due  to  the  presence  of  the  crack  can  be 
seen.  Similarly,  a  comparison  of  Figure  29  with  Figures  15  and  19  shows  the  effect  along  the 
lower  interface.  Both  top  and  bottom  interfaces  show  a  compressive  <jy  in  the  crack  flank 
region.  In  the  case  of  the  upper  adherend,  the  general  state  of  compression  in  the  adhesive  in 
the  crack  flank  may  be  due  to  the  relatively  stiff  adherend  pulling  the  otherwise  relaxed 
material  up  against  the  still  intact  adhesive  material.  In  the  lower  adherend  (Figure  29),  the 
compression  may  be  a  secondary  effect:  due  to  the  larger  longitudinal  stress  in  the  region  than 
in  the  uncracked  case,  the  discontinuity  in  material  properties  causes  a  slight  compression. 

It  is  significant  that  the  displacement  solution  shows  the  crack  is  actually  open.  This  is 
sensible   since   the   uncracked   solution   shows  a  tensile  stress  or   peeling   stress.     It    is   also 
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significant  to  note  that  the  stress  intensity  factors  are  of  the  same  order  of  magnitude,  with  k,, 
actually  being  somewhat  larger  than  k,  (see  Table  4).  This  indicates  the  failure  is  mixed 
mode,  rather  than  opening  mode  only,  as  cited  in  some  experimental  work  [25].  Again  this  is 
sensible,  physically. 

A  further  implication  from  current  results  may  indicate  that  if  the  adhesive  material  is 
sufficiently  ductile  to  accommodate  the  singularity  at  the  crack  tip  by  local  plastic  deformation 
rather  than  by  brittle  fracture,  then  the  interfacial  peel  stress  and  the  shear  stress  values  could 
be  lower  than  the  end  of  lap  stresses  for  the  uncracked  problem.  This  indicates  that  the 
perturbation  distance  of  the  stress  intensity  at  the  crack  tip  is  no  greater  than  the  perturbation 
distance  of  the  stress  concentration  due  to  the  geometric  discontinuity  of  the  uncracked 
problem,  which  then  implies  that  the  perturbation  distance  is  more  strongly  a  function  of 
material  properties  (as  seen  in  Section  V)  than  of  geometry. 

Section  VII:  Summary  and  Conclusions 

The  first  part  of  the  research  examined  the  solution  of  the  single  lap  shear  adhesive 
joint  with  no  crack.  The  improvement  of  the  results  with  respect  to  earlier  efforts  is  due  to  two 
factors.  First,  assumed  stress  hybrid  elements  were  used  in  most  of  the  critical  areas  of  the 
problem.  Since  the  hybrid  element  yields  the  same  degree  of  accuracy  as  a  second  order 
displacement  model  [21],  it  is  inherently  more  accurate  than  the  constant  strain  elements  used  in 
previous  numerical  solutions  [10,11]  or  "classical"  solutions  which  do  not  address  stress  in  the 
adhesive  at  all.  Second  is  the  use  of  more  elements  across  the  adhesive  layer  to  determine  the 
variation  in  stress  across  the  layer  and  to  obtain  a  good  approximation  of  the  value  of  the 
stress  at  the  interface. 
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Three  conclusions  can  be  drawn  from  the  first  part  of  the  research.  First,  there  is  a 
strong  stress  gradient  across  the  thickness  of  the  adhesive  as  well  as  along  the  lap,  particularly 
in  the  vicinity  of  the  end  of  the  lap.  Second,  no  solution  to  the  lap  shear  problem  can  ignore 
the  longitudinal  stress  in  the  adherend  as  is  commonly  done  in  the  literature  where  only  the 
peeling  stress  and  the  shear  stress  are  discussed.  The  importance  of  the  longitudinal  stress  has 
been  shown  experimentally,  but  it  has  never  before  been  specifically  shown  numerically  or 
analytically.  Third,  the  magnitude  of  the  stress  at  the  end  of  the  lap  is  higher  than  previously 
reported.  This  is  a  result  of  the  assumption  of  an  infinitessimal  adhesive  thickness  in  the 
analytical  solutions  and  due  to  the  use  of  constant  strain  elements  and  a  mesh  that  was  too 
coarse  in  earlier  numerical  solutions. 

The  second  part  of  the  research  is  the  examination  of  the  adhesive  bonded  single  lap 
joint  with  a  short  edge  crack.  An  assumed  stress  hybrid  super  element  was  used  at  the  crack 
tip  which  includes  terms  capable  of  describing  the  singular  behavior  at  the  crack  tip.  It  is 
significant  that  this  is  the  first  known  fully  operational  computer  program  utilizing  the  general 
nine  node  super  element  and,  therefore,  the  first  application  of  the  element  to  a  real  problem, 
and  the  first  analytical  attempt  to  address  the  adhesive  fracture  problem  with  real  physical 
dimensions  and  material  properties. 

The  results  from  the  second  part  were  not  thoroughly  investigated  and  are  preliminary 
in  nature.  They  indicate  that  the  stress  intensity  factors  are  related  to  the  modulus  ratio  by  a 
power  law.  They  also  show  that  the  stress  intensity  factor  for  inplane  shear  is  not  only  of  the 
same  order  of  magnitude  as  the  stress  intensity  factor  for  the  opening  mode,  but  actually 
slightly  larger.   The  comparability  in  the  magnitude  of  kj  and  kn  is  contraiy  to  the  conclusions 
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drawn  in  experimental  work  such  as  Mostovoy  [24]  and  others  which  indicates  kj  as  being  an 
order  of  magnitude  larger  than  kn.  The  fact  that  most  polymeric  adhesives  are  more  sensitive 
to  fracture  in  the  opening  mode  makes  experimental  determination  of  inplane  effects  difficult, 
which  would  account  for  experimental  underestimation  of  kn.  The  analysis  shows  that  stresses 
are  relieved  in  the  crack  flank  region  and  a  slight  compression  exists  normal  to  the  interface. 
Finally,  the  results  from  the  analysis  of  the  lap  joint  with  a  short  crack  indicate  that  there  is  no 
effect  on  the  uncracked  solution  in  the  "far"  field,  and  that  far  can  be  interpreted  as  being  of 
the  order  of  5-10  adhesive  thicknesses  along  the  lap  away  from  the  crack  tip. 

Section  VIII:  Suggestions  for  Further  Work 

The  results  of  the  research  reported  warrant  further  study  of  the  adhesive  bonded  lap 
shear  joint.  The  most  promising  approach  to  the  problem  is  a  parametric  study  of  the  effect 
of  the  many  variables  involved,  first  in  the  uncracked  case  and  then  using  that  as  a  basis  in 
the  cracked  case.  Because  no  rigorous  study  of  the  uncracked  case  exists,  a  catalogue  of  results 
for  the  uncracked  case  must  be  built  up,  so  the  cracked  case  can  be  compared  to  it,  to  see  the 
effect  of  the  crack. 

A  parametric  study  on  the  uncracked  adhesive  bonded  joints  should  include  not  only 
single  lap  shear  joints,  but  also  double  lap  shear  and  scarf  joints.  For  these  the  effect  of  the 
following  geometric  variables  should  be  investigated,  in  addition  to  the  effect  of  different 
adhesive  thickness  presented  herein: 

1.  Overlap  length  -  particularly  with  respect  to  the  effect  on  longitudinal  stress 

2.  Shape  of  adhesive  at  end  of  lap 


<H 


3.  Non-uniformity  of  adhesive  layer 

4.  Thickness  of  adherends 

5.  Shape  of  adherends  in  overlap  region,  such  as  the  recommendations  of  Cheery 
and  Harrison  [36] 

Future  work,  must  also  study  the  effect  of  material  variables  in  addition  to  the  moduli 

* 

ratio  studied  herein  such  as: 

1.  Dissimilar  adherend  materials 

2.  Anisotropic  materials 

3.  Viscoelasticity 

A.   Material  incompressibility 
5.   Plasticity 

All  variables  mentioned,  except  the  last  three,  may  be  studied  with  the  present  computer 
program,  but  modifications  of  the  formulation  are  required  to  handle  any  of  the  last  three.  A 
systematic  study  should  also  include  joint  loadings  other  than  tensile,  such  as  bending  or 
biaxial  stress. 

When  results  are  in  hand  for  the  uncracked  lap  joint  are  determined,  all  of  the  variables 
mentioned  for  the  uncracked  case  should  be  studied  with  respect  to  various  crack  variables, 
such  as: 

1.  Crack  length 

2.  Edge  cracks 

3.  Internal   cracks 


45 

4.   Orientation  of  crack 

The  current  computer  program  can  handle  all  the  above  cases  with  the  noted  exceptions 
(viscoelasticity,  material  incompresibility,  and  plasticity).  The  only  problem  in  continuing  the 
research  is  manpower  and  computation  expense. 
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Appendix  1 
Coland  and  Reissner  Solution  to  Single  Lap  Shear  Adhesive  Joint  [1] 


where 


xy 


4L 


0L  cosh(j5x/t) 

(I*  3«)     

2t  sinh(0L/2t) 


+    3(1-0 


rxy  -  shear  stress  at  point  x 


T  =  load  on  lap  joint,  per  unit  width 


L  ■  length  of  the  overlap 


t  -  thickness  of  each  adherend 


x  »  distance  from  the  center  of  lap 


8Gnt 
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(l  rixT 


1  ♦  2  V2     tanh 


2t\     2Et 


G8  =»  shear  modulus  of  adhesive 
ta  =  thickness  of  adhesive  layer 


\  - 1  -  vl 


v  »  Poisson's  ratio  of  adherends 
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Mesh  Arrangement  in  Vicinity  of  End  of  Lap  with  no  Crack 

FIGURE  10 
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Mesh  Arrangement  in  Vicinity  of  End  of  Lap  with  Short  Edge  Crack 

FIGURE    11 
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Figure  21 
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Figure  22 
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Table  1 


Comparison  to  Goland   Reissner  Solution 


Distance  from 
end  of  lap  (in) 

Pee 

1  Stress  <Ty 

She 

ar  Stress  T„ 

/ 

Current 

Solut ion 

Current 

So  lu t ion 

C-R 

n i dp  lane 

Pt  B 

G-R 

Midp  lane 

Pt  B 

8.0005 

0.839 

0.832 

1.137 

0.599 

0.547 

0.735 

8.0015 

0.794 

0.791 

0.825 

0.585 

0.504 

0.612 

8.003 

0.730 

0.695 

0.725 

0.564 

0.485 

0.571 

8.006 

0.610 

0.541 

0.560 

0.524 

0.462 

0.520 

0.013 

0.371 

0.286 

0.310 

0.442 

0.394 

0.448 

0.034 

-0.055 

-0.012 

0.000 

0.265 

0.258 

0.266 

0.100 

-0.084 

-0.116 

-0.111 

0.054 

0.076 

0.077 

8.200 

0.006 

0.043 

0.046 

0.005 

0.001 

-0.004 

The  following  parameters  were  used  in  above  solutions 


Adherend  -  Aluminum  (E  =  107,  v  =  0.33) 


Elastic  Adhesive  (E  =  500,000,  v  -  0.35) 


Adhesive  thickness  =  0.005  inch 


Adherend  thickness  =  0.05  inch 


Overlap  length  =  0.5  inch 


Loading  - 1  psi  tension  in  adherend 
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Table  2 
Comparison  of  Stress  Intensity  Factor  k,  for  Symmetric  Loading  with  Existing  Solutions  [31] 

Stress  Intensity  Factor,  kj 


Method  of 

Evaluating  kj 

Direct 

Crack  Opening 

PianTong.Luk 

[30] 

Mesh       D.O.F. 

1x2         24 

0.840 

0.904 

2x4         72 

0.822 

0.797 

4x4        128 

0.812 

0.784 

4x6        184 

0.791 

0.790 

Five  Node  Symmetric 

Super  Element  [31] 

2x3         22 

0.793 

— 

4x5         58 

0.794 

— 

Nine  Node  General 

• 

Super  Element 

(Present  Solution) 

6x8        130 

0.793* 

— 

"Exact"  Solution 

Bowie  [33] 

0.79' 

1 

*kn  «  1.7  x  10"5 
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Table  3 
STIFFNESS  MATRIX  FOR  NINE  NODE  GENERAL  SUPER  ELEMENT  (18  X  18) 

Note:    Matrix  is  symmetric  so  only  lower  triangle  is  shown.   Part  1  is  columns  1-6,  Part  2 
is  columns  7-12,  and  Part  3  is  columns  13-18. 


0.6685E+88 

e.o 

8.2349E-G1 
-0.5313E-01 
-8.1687E+88 

8.14S9E+88 
-8.9582E-81 

0.63G3E-81 

PART  1-G.9480E-81 

-0.5295E-02 

-8.948CE-81 

8.5295E-82 
-8.9582E-81 
-8.63C3E-01 
-8.1667E+00 
-8.1459E408 

0.2349E-G1 

8.5313E-81 


8.5966E+08 

8.8809E-01 
-8.1196E+08 

0.83C3E-81 
-8.1G16E+0G 

8.2482E-81 
-8.1779E-81 

8.1039E-01 

8.6874E-83 
-8.1839E-81 

8.G874E-03 
-8.2482E-81 
-8.1779E-81 
-0.8383E-O1 
-0.1616E+00 
-0.8809E-01 
-8.1196E+00 


0.3252E+G0 
8.8209E-01 
-8.1210E+00 
-8.7351E-01 
-8.8528E-81 
8.151GE-81 
-0.7393E-G1 
-0.1007E-01 
-0.1978E-81 
-8.G995E-02 
-8.2882E-81 
-8.1271E-81 
-8.1491E-81 
-0.6485E-81 
-8.1294E-81 
-8.1723E-01 


G.3249E+00 

8.7G91E-81 

8.1380E-01 

-8.1766E-81 

-8.3225E-81 

-0.4102E-G1 

0.1H1E-O1 

-8.2088E-G1 

8.2428E-82 

-8.2641E-01 

-8.2235E-81 

-8.1711E-81 

-8.9991E-81 

G.1723E-8i 

-0.7848E-01 


8.5710E+88 
-8.1165E-81 
-8.1024E+G8 
-8.1546E+88 
-8.16G9E+GB 
-8.2484E-81 
-8.2777E-02 

8.8386E-G2 
-8.5198E-G2 
-8.1539E-01 

8.4875E-82 

8.2836E-81 
-8.1491E-81 

8.1712E-61 


8.5G89E+08 

8.623GE-01 
-8.2353E-01 
-G.1S29E+0G 
-8.1451E+88 

8.4531E-B2 
-G.3251E-02 
-8.1938E-81 
-8.2785E-01 
-8.2036E-O1 
-0.1215E4O8 

0.6484E-81 
-0.9991E-81 


8.3047E+08 
-8.3814E-81 

8.2251E-01 

-8.7893E-01 

-8.1281E-81 

PART  2-8.7781E-82 

-8.5711E-82 

8.1552E-82 
-8.5199E-82 

8.1929E-01 
-8.2C82E-81 

8.2641E-81 


8.1321E+00 
8.7208E-01 

-8.1191E-01 
8.7961E-G2 
8.1884E-01 

-8.I552E-82 

-8.1528E-81 
8.1539E-01 

-8.2785E-01 
8.1271E-01 

-8.2235E-01 


8.358SE+00 
8.1875E+00 
-0.1G90E-81 
-0.4453E-82 
-0.1281E-O1 
-8.7961E-02 
-0.1777E-02 
-8.4530E-82 
-0.1976E-81 
8.2O88E-01 


8.1414E+00 

8.4453E-82 

-8.1444E-81 

8.7702E-82 

8.1884E-81 

-8.8307E-82 

-8.3252E-02 

8.6995E-02 

8.2429E-82 


8.3585E+08 
-8.1075E+88 

8.2251E-01 
-8.7288E-G1 
-8.16G9E+B8 

8.1529E+08 
-8.7393E-81 

8.4I02E-81 


8.1414E+00 
8.7S93E-81 

-0.1191E-81 
8.2483E-81 

-8.1451E+O0 
8.1887E-81 
8.1141E-81 


Table  3  Continued 


0.3O47E+OO 

0.3014E-O1 

PART  3-0.1O24E+0O 

-0.623GE-81 

-8.8528E-01 

0.1766E-01 


0.1321E+00 

0.1546E+OO 

-0.2353E-O1 

-0.1516E-01 

-8.3225E-01 


0.571OE+OO 

0.1165E-O1 

-0.121OE+GO 

-8.7G91E-01 


0.S689E+OO 
0.7351E-01 
0.138OE-O1 


0.3252E+OO 
-0.8209E-01 


0.3249E+OO 
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Table  4 


Stress  Intensity  Factors  kj  and  kn  for  Selected  Moduli  Ratios 


Adherend  to 

Stress  Intensity  Factor 

Adhesive  Modulus 

I 

ki 

kn 

.1040 

.1096 

2 

.0774 

.0845 

20 

.0253 

.0291 

200 

.0083 

.0101 

Geometry  as  shown  in  Figure  1 


Mesh  arrangement  as  shown  in  Figure  11 


Total  Decrees  of  Freedom  =  1144 


Total  Number  of  Elements  =  509 
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